#include "AnaDef.h"
#include "AnaAliDef.h"
#include "TTreeStream.h"
#include "AliTRDdEdxBaseUtils.h"
#include "AliTRDdEdxCalibUtils.h"
#include "AliTRDdEdxReconUtils.h"
#include "AliTRDtrackV1.h"
#include "AliTPCseed.h"
#include "IniRun.h"
#include "AliTRDCalTrapConfig.h"
#include "AliTRDcalibDB.h"
#include "AliESDtrack.h"

TTree *ftt = 0x0;

Int_t fopt=-999;
Int_t ftbmin = -1;
Int_t ftbmax = 1000;

TString fcdbpath;

#define EPSILON 1e-12

void analysis()
{
  Int_t run = -999;
  AliTRDtrackV1 *trdtrack=0x0;
  AliTPCseed *tpcseed =0x0;
  AliESDtrack * esdtrack = 0x0;
  Int_t sector=-999, istack=-999;
  Double_t ltbg=-999;

  Double_t tpcsig=-999;
  Int_t tpcncls = -999;
  Int_t pid=-999;

  Int_t trig=-999, itrk=-999, mag=-999, charge=-999;

  Double_t cen = -999;

  ftt->SetBranchAddress("run", &run);
  ftt->SetBranchAddress("trdv1", &trdtrack);
  ftt->SetBranchAddress("isector", &sector);
  ftt->SetBranchAddress("istack", &istack);
  ftt->SetBranchAddress("tpcsig", &tpcsig);
  ftt->SetBranchAddress("tpcseed", &tpcseed);
  ftt->SetBranchAddress("esdtrack", &esdtrack);
  ftt->SetBranchAddress("mag", &mag);
  ftt->SetBranchAddress("pid", &pid);
  ftt->SetBranchAddress("charge", &charge);
  ftt->SetBranchAddress("ltbg", &ltbg);

  ftt->GetEntry(0);

  //=============================

  const TString rtype = AliTRDdEdxBaseUtils::GetRunType(run);
  if(rtype.Contains("COSMIC")){
    ftt->SetBranchAddress("itrk", &itrk);
    ftt->SetBranchAddress("trig", &trig);
  }
  else if(rtype.Contains("PP")){
  }
  else if(rtype.Contains("PBPB")){
    ftt->SetBranchAddress("pbcen",&cen);
  }
  else{
    printf("newdraw.C run type error!! %s %d\n", rtype.Data(), run); exit(1);
  }

  const TString ryear = rtype(0,4);
  //void SetOCDB(const Int_t kgrid, const Int_t kraw, const Int_t year, const Int_t xxrun, const Bool_t k1kg)
  AliCDBManager *cdb = SetOCDB(0, 0, ryear.Atoi(), run, rtype.Contains("1KG"));
  if(fopt==11){
    cdb->SetSpecificStorage("TRD/Calib/PHQ",Form("local://%s", fcdbpath.Data()));
  }
  cdb->Print();

  //=============================================================
  //=============================================================
  TList *list = 0x0;
  TTreeSRedirector *outStream = 0x0;
  if(fopt==10){
    list = new TList;
    list->SetOwner(kFALSE);
    AliTRDdEdxCalibUtils::IniHistArray(list,kFALSE);
  }
  else if(fopt==11){
    AliTRDdEdxCalibUtils::SetObjArray(AliTRDcalibDB::Instance()->GetPHQ());
  }

  //=============================================================
  //=============================================================

  const Int_t nevt = ftt->GetEntries();
  for(Int_t ievt=0; ievt<nevt; ievt++){
    if(ievt%10000==0) printf("Entry: %d/%d\n", ievt, nevt);

    ftt->GetEntry(ievt);

    //can not use tpcncls<0 to judge! since it is not -999 once set by tpcseed
    if(tpcseed){
      tpcncls = tpcseed->GetNumberOfClusters();
    }

    Int_t itype = -999;
    if(rtype.Contains("COSMIC")){//cosmic
      itype = itrk;
    }
    else if(rtype.Contains("PP")){//pp
      itype = 2;
    }
    else{//PbPb
      itype = 10 + Int_t(cen/10.);
    }

    if(itype==0 || itype==1){
      /*
      //only for cosmic ---->
      const TString sidesec(" 0 8 9 17 ");
      if(sidesec.Contains(Form(" %d ", sector))){
        continue;
      }
      //<---
      */
    }
    
    //------------------------------------------
    //------------------------------------------

    if(fopt==10){
      Double_t scale = 1;

      if(itype==0 || itype==1){
        //cosmic
        if(ltbg<1.65 || ltbg>2.5) //for feb up to 1.6 still see electron
          continue;
      }
      else{
        scale = AliTRDdEdxCalibUtils::GetCalibTPCscale(tpcncls, tpcsig);

        if(scale<0)
          continue;

        // if(ltbg<0.5 || ltbg>2)
        //   continue;
      }

      if(AliTRDdEdxBaseUtils::IsExBOn()){
        AliTRDdEdxCalibUtils::FillHist(trdtrack, 0, mag, charge, scale);
        AliTRDdEdxCalibUtils::FillHist(trdtrack, 1, mag, charge, scale);
      }
      else{
        AliTRDdEdxCalibUtils::FillHist(trdtrack, 0, -1, -1, scale);
        AliTRDdEdxCalibUtils::FillHist(trdtrack, 1, -1, -1, scale);
      }
    }
    else if(fopt==11){
      if(itype>=2){
        //only for pp--->
        //cut away overlap of pion and proton in pp
        if(pid==4 && ltbg>0.4)
          continue;
        if(pid==2 && ltbg<0.5)
          continue;
        
        if(pid!=0 && pid!=2 && pid!=4)
          continue;
        //<---
      }

      TVectorD trdarrayQ0(200), trdarrayX0(200), trdarrayQ1(200), trdarrayX1(200);
      Int_t trdncls0 = -999, trdncls1 = -999;
      Int_t trdnch0 = -999, trdnch1 = -999;
      Double_t trdQ0 = -999, trdQ1 = -999;

      //CookTruncatedMean(0, mag, charge, kTRUE, fNchamberdEdx, fNclusterdEdx);

      if( AliTRDdEdxBaseUtils::IsExBOn() ){
        trdQ0 = trdtrack->CookTruncatedMean(0, mag, charge, kTRUE, trdnch0, trdncls0, &trdarrayQ0, &trdarrayX0, ftbmin, ftbmax);
        trdQ1 = trdtrack->CookTruncatedMean(1, mag, charge, kTRUE, trdnch1, trdncls1, &trdarrayQ1, &trdarrayX1, ftbmin, ftbmax);
      }
      else{
        trdQ0 = trdtrack->CookTruncatedMean(0, -1, -1, kTRUE, trdnch0, trdncls0, &trdarrayQ0, &trdarrayX0, ftbmin, ftbmax);
        trdQ1 = trdtrack->CookTruncatedMean(1, -1, -1, kTRUE, trdnch1, trdncls1, &trdarrayQ1, &trdarrayX1, ftbmin, ftbmax);
      }
      //------
      if((itype==0 || itype==1) && (ltbg>1.65 && ltbg<2.5)){
        //cosmic
        static Int_t kprint = 100;
        if(kprint<0){
          printf("\nsummary:\n");
          
          printf("trd0 %d %15f\n", trdncls0, trdQ0); for(Int_t iq=0; iq<trdncls0; iq++){printf("trdarrayX0[%3d] = %15.0f; trdarrayQ0[%3d] =  %15f;\n", iq, trdarrayX0[iq], iq, trdarrayQ0[iq]);}
          printf("trd1 %d %15f\n", trdncls1, trdQ1); for(Int_t iq=0; iq<trdncls1; iq++){printf("trdarrayX1[%3d] = %15.0f; trdarrayQ1[%3d] =  %15f;\n", iq, trdarrayX1[iq], iq, trdarrayQ1[iq]);}
          printf("\n");
        }
        kprint++;
      }

      //----------------------

      Int_t rawnch  = trdtrack->GetNumberOfTracklets();      
      Int_t rawncls = trdtrack->GetNumberOfClusters();
      Double_t pTRD = TMath::Power(10, ltbg)*AliPID::ParticleMass(pid);

      Double_t eta = esdtrack->Eta();
      Double_t phi = esdtrack->Phi();

      const AliTRDseedV1 * seed0 = AliTRDdEdxBaseUtils::GetFirstTracklet(trdtrack);
      const AliTRDseedV1 * seed1 = AliTRDdEdxBaseUtils::GetLastTracklet(trdtrack);

      Double_t dydx0 = AliTRDdEdxBaseUtils::Getdydx(seed0);
      Double_t dydx1 = AliTRDdEdxBaseUtils::Getdydx(seed1);

      Double_t dzdx0 = AliTRDdEdxBaseUtils::Getdzdx(seed0);
      Double_t dzdx1 = AliTRDdEdxBaseUtils::Getdzdx(seed1);

      if(!outStream){
        outStream = new TTreeSRedirector(Form("outpar%d_%010d.root", fopt, run));
      }
      
      (*outStream)<<"newdraw_Tree"<<
        "run="   <<run<<
        "itype="  <<itype<<
        "trig="  <<trig<<
        "cen="<<cen<<
        "mag="   <<mag<<
        "charge="<<charge<<
        "ltbg="  <<ltbg<<
        "pid="   <<pid<<
        "mom="   <<pTRD<<
        "isector="<<sector<<
        "istack="<<istack<<
        
        "ievt="  <<ievt<<
        
        "rawnch="<<rawnch<<
        "rawncls="<<rawncls<<
        
        "trdnch0=" <<trdnch0<<
        "trdncls0=" <<trdncls0<<
        "trdQ0=" <<trdQ0<<

        "trdnch1=" <<trdnch1<<
        "trdncls1=" <<trdncls1<<
        "trdQ1=" <<trdQ1<<
        
        "dydx0="<<dydx0<<
        "dydx1="<<dydx1<<

        "dzdx0="<<dzdx0<<
        "dzdx1="<<dzdx1<<

        "eta="<<eta<<
        "phi="<<phi<<

        "tpcsig="<<tpcsig<<
        "tpcncls="<<tpcncls;
      
      if(rtype.Contains("COSMIC")){
        (*outStream)<<"newdraw_Tree"<<
          "trdarrayQ0="<<&trdarrayQ0<<
          "trdarrayX0="<<&trdarrayX0;
      }
      
      (*outStream)<<"newdraw_Tree"<<
        "\n";
    }
  }
  
  if(fopt==10){
    //list->ls();
    printf("Start CalibOutput...\n");

    TList *lout = new TList;
    lout->SetOwner();

    TTreeSRedirector *calibStream = new TTreeSRedirector(Form("TRDCalibStream_%010d.root", run));

    for(Int_t iter=0; iter<AliTRDdEdxCalibUtils::GetHistArray()->GetSize(); iter++){
      THnBase *hi = (THnBase*) AliTRDdEdxCalibUtils::GetHistAt(iter);
      TObjArray *obji = AliTRDdEdxCalibUtils::HistToObj(hi, run, lout, calibStream);
      //printf("test analyze %s\n", obji->GetName());
      AliTRDdEdxCalibUtils::GetObjArray()->AddAt(obji, iter);
    }
    
    TFile *fout=new TFile(Form("TRDCalibList_%010d.root", run),"recreate");
    fout->cd();
    list->Write();
    lout->Write();
    fout->Save();
    fout->Close();
    delete fout;

    delete calibStream;
    delete lout;
    
    AliTRDdEdxCalibUtils::GenerateOCDB(run);

    printf("CalibOutput Ends!\n");    
    //list->ls();

    delete list;
    AliTRDdEdxCalibUtils::DeleteHistArray();
  }
  else if(fopt==11){
    AliTRDdEdxCalibUtils::DeleteObjArray();
    delete outStream;  
  }
  
}

void newdraw(const TString fin, const Int_t opt, const Int_t tb0 =-1, const Int_t tb1 =10000, const TString cpath="")
{
  printf("fin: %s\n", fin.Data());

  fopt = opt;

  ftbmin = tb0;
  ftbmax = tb1;
  
  fcdbpath = cpath;

  printf("\n*********** opt %d, tb %d %d **************\n\n", fopt, ftbmin, ftbmax);
  //=====================
  //=====================

  AliTRDdEdxBaseUtils::PrintControl();

  TStopwatch timer;
  timer.Start();

  TFile *ftmp = new TFile(fin);
  ftt=(TTree*)gDirectory->Get("tree");
  if(!ftt){
    printf("no tree!!\n");
    gDirectory->ls();
    exit(1);
  }
  
  analysis();

  timer.Stop();
  timer.Print();

  delete ftmp;
}

int main(const int argc, char *argv[])
{
  int nums[argc];
  for(int ii=0; ii<argc; ii++){
    printf("%d: %s\n", ii, argv[ii]);
    nums[ii] = atoi(argv[ii]);
  }

  if(argc!=6){
    printf("argc!!! %d\n", argc);
    return 1;
  }

  newdraw(argv[1], nums[2], nums[3], nums[4], argv[5]);
  return 0;
}
